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Abstract. The excitation of internal gravity waves by penetrative convective plumes is investigated using 2-D direct 
simulations of compressible convection. The wave generation is quantitatively studied from the linear response 
of the radiative zone to the plumes penetration, using projections onto the g-modes solutions of the associated 
linear eigenvalue problem for the perturbations. This allows an accurate determination of both the spectrum and 
amplitudes of the stochastically excited modes. Using time-frequency diagrams of the mode amplitudes, we then 
show that the lifetime of a mode is around twice its period and that during times of significant excitation up to 
40% of the total kinetic energy may be contained into g-modes. 
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1. Introduction 

Although their detection in the spectrum of solar oscilla- 
tions has been not clearly confirmed (e.g., Turck-Chieze 
et al. 2004; Gabriel et al. 2002), internal gravity waves 
(hereafter IGWs) propagating in the radiative zones of 
late-type stars have recently been invoked in attempts to 
explain the Li abundance of cool stars and the rigid rota- 
tion of their radiative interiors. 

The former problem is also referred to as the Li- 
dip problem and concerns the dependence of the lithium 
abundance on the spectral type for some main-sequence 
stars. Models based on the extension of the surface con- 
vection zone down to the nuclear burning region (Ibcn 
1965) and models which take into account the transport 
of Li both by meridional circulation and shear-induced 
turbulence (Talon & Charbonnel 1998) are indeed not 
quite satisfactory to reproduce the Li-dip. The latter prob- 
lem concerning the rigidity of the Sun's radiative inte- 
rior is most clearly revealed by helioseismology (Brown 
et al. 1989). Both rotation-induced turbulent diffusion 
and wind-driven meridional circulation fail to redistribute 
enough angular momentum over the lifetime of the Sun to 
rotate rigidly (Zahn et al. 1997). Likewise, the hypothe- 
sis of a large-scale poloidal magnetic field leads to prob- 
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lems, because it may transmit under certain circumstances 
the differential rotation of the convection zone to the core 
(owing to Ferraro's law; e.g., MacGregor & Charbonneau 
1999). 

Zahn (1994) showed that the two problems are coupled 
and that they should to be explained by a single model. 
Being inspired by meteorological studies of the wave trans- 
port taking place in the Earth stratosphere (Bretherton 
1969; Alexander & Pfister 1995), Schatzman (1993) was 
the first to propose internal gravity waves as an efficient 
transport mechanism in stellar radiative zones. Later, this 
idea was tested by many authors (Zahn et al. 1997; Kumar 
& Quataert 1997; Kumar et al. 1999; Talon, Kumar & 
Zahn 2002) who showed that internal gravity waves trans- 
port momentum on a rather short timescale such that the 
rotation of the solar core becomes nearly uniform. A re- 
maining problem is the excitation of these deep gravity 
waves since, unlike pulsating white dwarfs, a K-mechanism 
based on hydrogen and helium ionization zones is not ap- 
plicable here. 

The most wide-spread excitation model is based on 
penetrative convection from neighboring convection zones. 
Strong downward plumes are known to extend a substan- 
tial distance into the adjacent stable zones so that internal 
gravity waves can be randomly generated. 2-D and 3-D 
direct numerical simulations of superposed stable and un- 
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stable layers have confirmed this scenario since gravity-like 
waves have been well observed in stable regions (Hurlburt 
et al. 1986, 1994; Brandenburg et al. 1996; Brummell et 
al. 2002). 

The aim of this paper is to investigate in detail the ex- 
citation of IGWs by overshooting in an high-resolution 
2-D hydrodynamical simulation of a three-layer model, 
consisting of a convective zone (hereafter CZ) embed- 
ded between an upper cooling zone and a lower stably 
stratified radiative zone (hereafter RZ). In all previous 
numerical studies, gravity waves have been detected us- 
ing Fourier's analysis such as, for example, the (k, oj)- 
diagram. However, the link with the eigenmodes of sta- 
ble regions was not investigated. In particular, it has not 
clearly been demonstrated that the excited waves really 
correspond to gravity waves since the only criterion was 
the Brunt- Vaisala, frequency as an upper bound. 

In this work, g-modes are rigorously measured using 
the method that we have presented and tested on the 
gravity mode oscillations of an isothermal atmosphere in 
Dintrans & Brandenburg (2004, hereafter Paper I). Our 
method mainly relies on the projection of the velocity held 
of the simulation onto the basis shaped from the solutions 
of the associated linear eigenvalue problem for the pertur- 
bations; i.e., the theoretical g-modes of the RZ. Hence the 
mode amplitudes are simply given by the time-dependent 
basis coefficients, which allows a quantitative study of the 
excitation mechanism. In other words, we investigate the 
generation of IGWs by penetrative convection from the 
linear response of the RZ to this penetration. 

We begin by presenting our hydrodynamical 2-D model 
consisting in three superposed polytropic layers, and give 
some details on the code we use to solve it numerically 
(Sect. We then give the main properties of the ob- 
tained simulation of penetrative convection and show that 
the classical detection method of gravity mode oscillations 
based on Fourier's transforms in both space and time fails 
to give reliable results on this problem (Sect.[3J). We then 
introduce our detection technique based on the anelastic 
subspace and apply it to find the properties of the g-modes 
propagating in the numerical simulation (Sect.^J. Finally, 
we conclude in Sect. El by giving some outlooks of this 
work. 

2. The hydrodynamical model 

2.1. The basic equations 

We adopt Cartesian coordinates (x, z) where x denotes the 
horizontal direction and z is depth pointing downward as 
the gravity g. Our system is composed of a convection 
zone of depth d — z% — Z2, embedded between two stable 
layers (Figure QJ. We assume that the gas is monatomic 
and perfect, so its equation of state is given by 

p = (7 — l)pe with 7 = c p /c v = 5/3, 

where p is the pressure, p the density, e the internal energy, 
and 7 is the ratio of specific heats c p and c v . 
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Fig. 1. Geometry of the computational domain. 



We solve the following set of hydrodynamical equations 
(conservation of mass, momentum and energy): 



f Dlnp 
Dt 



divit, 



= -(7 - 1) (Ve + eVlnp) + g + - V • {puS), (1) 
Dt p 

De 1 

_ = _( 7 _i)edivu+ -V ■ UCVe) + 2uS 2 + Q, 
Dt p 



where u is the velocity and D/Dt = d/dt + u ■ V is the 
total derivative. In addition, v denotes the constant kine- 
matic viscosity 1 and K. = K/c v the radiative conductivity 
divided by c v 2 . To reproduce the radiative cooling taking 
place at a star's surface, we add a cooling term Q in the 
energy equation given by 



Q 



(e-etop)r x (z), 



(2) 



where r _1 (z) is the cooling rate profile which is set 
equal to zero everywhere except close to the surface, i.e. 
r _1 (z) 7^ only for z = [zi,z 2 ]- Because of this efficient 
cooling, the surface layer tends to be isothermal and the 
internal energy e is almost constant there. Finally, S de- 
notes the traceless strain tensor given by 



duj 
dx 3 



duj 
dx,. 



- -SijV-u. 



1 In Hurlburt et al. (1986, 1994), Cattaneo et al. (1991) or 
Bogdan et al. (1993), it is the dynamical viscosity /i = pv 
which is constant. However, the kinematic viscosity becomes 
very large as the density tends to zero leading to a highly 
viscous surface layer. 

2 Hereafter, we will simply refer to K. as the radiative con- 
ductivity. 
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2.2. Boundary conditions and the initial setup 

We assume that In p, u and e are periodic in the horizontal 
direction and adopt the following conditions at the upper 
and lower boundaries: 

du dw 



e = etop at z = z\ and — fixed at z = 24, 
\ oz 

where u and w correspond to the horizontal and vertical 
velocities, respectively. 

Following Brandenburg et al. (1996), we choose the 
depth of the unstable layer d as the unit of length, (d/g) 1 / 2 
as the unit of time and the initial value pbot of the density 
at the bottom of the convection zone (hereafter BCZ) as 
the unit of density [velocities are thus measured in units of 
\fgd, i.e. the free-fall velocity of the unstable layer divided 
by y/2, and fluxes in units of p\ Jot (gd) 3 ^ 2 ]. Finally, the 
dimensionless gravity is set equal to unity in the whole 
domain, i.e. 5 = 1 everywhere. 

The initial state is computed using polytropes in hy- 
drostatic and radiative equilibrium (e.g., Hurlburt et al. 
1986). In Appendix 1X1 we give the details of the initial 
setup, in particular the mixing length solution in the cal- 
culation of the CZ stratification, which helps to accelerate 
the numerical convergence towards the thermally relaxed 
state. 

2.3. The numerical method 

We use the hydrodynamical code described in Nordlund 
& Stein (1990) to advance the fully nonlinear set of Eqs. 
Q. In this code, spatial derivatives are computed using 
sixth order compact finite differences (Lele 1992) whereas 
the time advance is performed using a third order explicit 
Hyman scheme (Hyman 1979). Corresponding timesteps 
are usually 25% of the Courant-Friedrich-Levy timestep 
defined by 

StcFh = Az/ max(c s , \u\, av/Az, ax/Az), 

where c s is the sound speed, \ — Kji^p) the radiative 
diffusivity and a ~ 4 an empirical factor determining the 
length of the diffusive timestep. In practice, it is the radia- 
tive diffusion term which limits the timestep due to both 
the high vertical resolution and large radiative diffusivities 
of radiative zones. 

Contrary to acoustic modes, gravity modes correspond 
to waves with long periods, of about twenty time units in 
our numerical simulations, typically. Hence, once achieved 
the thermal relaxation of the radiative zone, we still need 
to integrate the dynamical equations over very long times 
to capture so much periods as possible (one typically needs 
runs as long as one thousand time units). We then chose 
to concentrate on a particular 2-D simulation with an high 
spatial resolution 256 x 512 (i.e. 256 mesh points in the 
horizontal and 512 points in the vertical) and an aspect 
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Fig. 2. Snapshots of the velocity field superimposed on 
a grey scale representation of the internal energy fluc- 
tuations (i.e. or temperature) at three different times 
t = [900, 903, 906]. The horizontal dashed lines delimit 
the convection zone z = [0, 1] located above the radiative 
zone z — [1,3]. 

ratio A = L x /d = 4 (L x being the horizontal extent of 
the computational domain). The main properties of this 
simulation are the following: 

mi = —0.9, TO2 = —0.8, 7713 = 3, e top = 0.3, 
< Ftot = 5 x 10- 3 , v = 5 x 10- 3 , (4) 

Pr 2 = 12.5, Pr 3 = 0.625, 

where Pr = vj\ denotes the Prandtl number. The corre- 
sponding Rayleigh number of the unstable layer is Ra ~ 
8.5 x 10 5 , using the definition in Gough et al. (1976). The 
kinematic viscosity v cannot be too small for a given reso- 
lution as it should satisfy d/ Az ~ Re 3//4 , where Az is the 
smallest mesh interval in the vertical direction and Re is 
the Reynolds number (e.g., Landau & Lifshitz 1980). The 
value we chose is a reasonably "safe" one as this is around 
fifty times larger than the minimum viscosity f m in ~ 10 -4 
imposed by both the resolution 256 x 512 and the mean 
Reynolds number of our simulation. 
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Fig. 3. Normalized vertical profiles of the radiative (solid 
dark line), convective (dashed line) and kinetic (dot- 
dashed line) fluxes, with their sum (solid grey line). The 
vertical dotted lines denote the CZ limits. 
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Fig. 4. Evolution with time of the penetration extent A, 
with its time-averaged value (A) ~ 0.56 denoted by the 
horizontal dotted line. 

3. Nature of the penetrative convection and the 
IGW detection problem 

Since the pioneering numerical simulations of Hurlburt et 
al. (1986), the general features of penetrative compressible 
convection are well known and we are simply giving here 
the main properties of such a flow. 

Figures and [3] represent typical asymmetrical pat- 
terns and mean vertical radiative, convective and kinetic 
fluxes that we obtain in our numerical simulation of com- 
pressible convection with penetration, the fluxes being 
computed from their usual definition (e.g., Hurlburt et 
al. 1986) 



rad 



dz 



< F conv = -Cp(pwT') t , 



(5) 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 
depth 

Fig. 5. Mean vertical profile of the Peclet number, where 
dotted lines mark the CZ limits. 

where the overbars and brackets denote the horizontal and 
temporal means, respectively, and T" is the temperature 
fluctuation about the horizontal mean. Compressibility ef- 
fects destroy the usual symmetric velocity field observed in 
Boussinesq convection as downdrafts become stronger and 
more concentrated than the broader upward flow (Graham 
1975; Hurlburt et al. 1984). The main consequence is that 
some localized downward-directed narrow plumes appear 
which transport a significant kinetic flux in this direction 
(around 30% of the total flux at the BCZ in this simula- 
tion; see the dot-dashed line in Fig. EJt- 

The presence of a lower stable layer below the convec- 
tion zone allows these downward plumes to penetrate some 
distance into the RZ. This convective penetration has been 
theoretically investigated in the astrophysical context by 
Zahn (1991) who showed that it strongly depends on the 
value of the local Peclet number Pe = wL/x {w and L be- 
ing the typical vertical velocity and size of those motions, 
respectively). 

Following Brummell et al. (2002), we define the pen- 
etration extent A as the vertical distance from the BCZ 
where the horizontally-averaged kinetic flux decreases to 
1% of its maximum value and Fig. 0] shows an example 
of its evolution with time. The averaged penetration is 
(A) ~ 0.56, that is, of order the half-size of the CZ, which 
is typical of hydrodynamical simulations at low Peclet 
numbers. Indeed, Pe falls down very rapidly just below 
the CZ due to the large radiative diffusivity of the radia- 
tive zone. In polytropic models, the radiative diffusivity is 
proportional to 1 + m, where m is the polytropic index, 
see Eq. (|A.3|) . so a convective blob experiences a jump in 
the Peclet number when it crosses the interface between 
the CZ and RZ zones 



Pes 
Pe 2 



X2 
X3 



1 + m 2 
1 + to 3 



(6) 



F kin = --(wp(u 2 + W 2 )} t , 



Figure [5] shows the vertical profile of the Peclet number. 
The estimate of the Peclet jump across the CZ-RZ inter- 
face using Eq.© gives Pe 3 /Pe 2 ~ 0.05 or Pe 3 ~ Pe 2 /20, 
which is effectively the one in Fig. El where Pe3/Pe2 — 
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Fig. 6. Left: Gray scale representations in the [z, o;)-plane 
of the low-frequency part of the temporal power spectra 
for two degrees £ = [1,2]. Right: the resulting spectra after 
integrating over depth. 

13/250 = 0.05. As a result, the cold blob thermalizes 
very rapidly in the radiative zone and loses its identity 
compared to its environment: the buoyancy braking thus 
disappears and allows the blob to continue by inertia in 
the stable zone, leading to the observed large penetra- 
tion. (In solar-type stars the relevant extent is actually 
much smaller, because here the thermal conductivity is 
still too large and therefore the overall fluxes, and hence 
the amount of flux is too large.) 

3.1. Finding the g-modes: the problem of the random 
excitation 

As discussed in Paper I, the classical technique to detect 
the g-modes propagating in an hydrodynamical simula- 
tion consists first of taking horizontal Fourier transforms 
of the vertical mass flux pw{x, z, t) for every time step, 
to get pw(i, z,t) 3 . Second, one computes power spectra 
for the individual time series pw(£, z, t), to get pw{£, z, li), 
and plots the resulting power spectra at a given degree £ 
in a (z, w)-plane to highlight the "shark fin" peaks corre- 
sponding to the eigenmodes (e.g., Fig. 5 in Paper I). 

Figure El shows the result of applying this method to 
the detection of £ = 1, I = 2 and I = 3 g-modes propa- 
gating in the simulation. Low-frequency peaks appear in 
the RZ (1 < z < 3), mainly in the 1=1 diagram, but 
these peaks are not as well defined as in Paper I. Indeed, 
we have focused in the previous paper on the simpler case 
of the g-modes of an isothermal atmosphere where IGWs 
were excited by the vertical free oscillations of a cold bub- 



3 Here and in the following, we define the horizontal 
wavenumber as k x = (2ir/L x ) £ = (n/2) £, where I denotes the 
mode degree and is a non-zero positive integer, £ = [1, 2, . . .], 
as purely radial g-modes cannot exist (e.g., Turner 1973). 



ble: once emitted, each global mode "stayed in the box" 
during a time that depends on the efficiency of the dis- 
sipation, which was mainly due to a weak viscosity. As a 
consequence, large-scale nonradial g-modes survived over 
long times and were very well visible in the temporal spec- 
tra, both in the (z, o;)-plane and in its depth-integrated 
representation (see Fig. 3 in Paper I). 

The situation is clearly different with penetrative con- 
vection as we have now to deal with a random excitation of 
IGWs. Indeed, strong downward plumes are not station- 
ary structures: they are born in the upper layers of the 
convection zone, are accelerated by the Archimedes force 
during their CZ crossing and, finally, end their life in the 
overshoot region where they transfer a large amount of 
their stored kinetic energy to the stably stratified medium, 
resulting in an internal gravity wave field. As the birth of 
these plumes is random, the resulting forcing of the RZ 
wave field is itself a random one, exactly as an hammer 
which randomly strikes the upper part of the RZ. We then 
understand why the detection of these IGWs is more diffi- 
cult than the simpler case of free bubble oscillations: once 
emitted by a penetrating plume, a global mode first prop- 
agates in the radiative zone while being subjected to the 
viscous and radiative dissipations. However, if a second 
plume arrives shorter afterwards, the mode pattern may 
be destroyed, resulting in partial interference and the cor- 
responding frequency peak may disappear from the spec- 
trum. In other words, the lifetime of a mode is not simply 
related to the diffusive processes but also to the frequency 
of plume penetration (as we will show in Sect. rO|l . 

With this in mind, it is also clear that the classical 
detection method based on successive Fourier transforms 
both in space and time of the vertical mass flux is not 
well adapted to detecting IGWs in simulations with pen- 
etrative convection. The temporal Fourier transforms are 
computed over the whole simulation while IGWs are only 
present during a short time interval. This results in a mix- 
ing of wave events with non-wavy turbulent events such 
that the spectra lack a well defined frequency (Fig. [fjl 
right). 



4. Results 

4.1. IGW identification using the anelastic subspace 

We have presented in Paper I a new detection method of 
IGWs in hydrodynamical simulations which is based on 
the anelastic subspace and we will give here only its main 
stages. The idea is to project the simulated velocity field 
onto the basis built with the anelastic eigenvectors which 
are solutions of the following linear problem for the adi- 
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Fig. 7. Mean vertical profile of the Brunt- Vaisala fre- 
quency N in the RZ, the dotted line denoting the BCZ. 

abatic perturbations (see also Dintrans & Rieutord 2001; 
Rieutord & Dintrans 2002) 

dz dz 
k C = for z = zi, z 4 , 

where £ = (£ X i6z) denotes the Lagrangian displacement 
vector (the velocity field being u = iu)£), po the mean 
density (i.e., po = (p)t) and A'' 2 the square of the Brunt- 
Vaisala frequency given by 

fl d\ne f 1\ d\np ~ 

7 dz \ 7 / dz 

To derive the anelastic subset {7J), we first filtered out the 
acoustic waves in the governing equations for the linear 
perturbations, second assumed that the time dependence 
of normal modes is of the form exp(iwt) and, finally, that 
their horizontal dependence follows from the imposed pe- 
riodicity in this direction as 

£ x oc cos(k x x) and £ z oc sin^a;). 

The coupled differential equations Q form a generalized 
eigenvalue problem of the form 

M.Ai>tn = U£ n MB*l>£n, (8) 

where ipe n = (£, x ,£,z) T is the eigenvector of degree £ and 
radial order n associated with the eigenvalue u>f ' , while 
M.A and M.b denote two differential operators. 

The profile of the Brunt- Vaisala frequency N in the 
radiative zone is shown in Fig.0 The g- modes spectrum is 
bounded by its maximum value, max(A^) ~ 0.46, while the 
typical frequency of a low-degree and low-order g-mode is 
given by u>j> n ~ N ~ 0.35, hence a nondimensional period 
T ln = 2tt/^„ ~ 18 (e.g., Turner 1973). 
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Fig. 8. First three anelastic eigenvectors if)g n — (t; x ,£,z) T 
at I = 1 and radial orders n = (solid line), n = 1 (dot- 
ted line) and n = 2 (dashed-dotted line) . Upper and lower 
panels show the vertical and horizontal displacements, re- 
spectively. Eigenvectors have been normalized by impos- 
ing fx* Polfanl 2 dz = 1. 

Figure [S] shows the vertical profiles of the first three 
anelastic eigenvectors of degree I = 1 and radial orders 
n = [0, 1, 2]. The associated (dimensionless) eigenvalues 
arewio = 0.278, wn = 0.184 and W12 = 0.133. Asg-modes 
are evanescent in a convectively unstable layer, each eigen- 
vector is trapped in the bottom stably stratified zone and 
its amplitude rapidly decreases in the convection zone, as 
observed for z < 1 where both £ x and £ z stop oscillating 
and tend to zero. 

Once the anelastic eigenvectors xpi n are computed, we 
can determine the amplitudes of the g-modes propagating 
in our simulation. Indeed, we have showed in Paper I that 
these eigenvectors are orthogonal to each other and form 
an orthogonal basis onto which the simulated velocity field 
can be projected as 

oo 

ue(z,t) = y^ J {il>tn,ut)il>tn + "rest", (9) 

71=0 

where the symbol ( , ) denotes a scalar product defined 

by 

(lpl,4>2) = V>1 • i>2 P0 dz, 

and the term "rest" contains all velocity components that 
are not due to IGWs (e.g., acoustic waves, convective ve- 
locities, etc). Here iig denotes the velocity field of degree 
£, that is, the horizontal Fourier transform of the veloc- 
ity field u(x, z, t). We then define the amplitude of the 
g-mode of degree I and order n as the time-dependent 
complex coefficient 

Cln(t) = (lp£n,U e ) 7 (10) 
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Fig. 9. Evolution with time of the real (solid lines) and 
imaginary (dotted lines) parts of the amplitude ci n of the 
g-modes of degrees I = [1, 2] and orders n = [0, 1]. 
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Fig. 10. Corresponding temporal power spectrum of the 
mode amplitude Q n (t) used in Fig. El f° r ^ = 1 and n = 0. 
The vertical dotted line marks the location of the theoret- 
ical frequency lu\q of the underlying g-mode. Amplitudes 
have been magnified by 10 5 for clarity. 

which corresponds to the basis coefficient for the anelastic 
eigenvector ipi n . 

4.2. The evolution of the mode amplitudes 

We show in Fig. [5] the resulting amplitudes obtained by 
applying our projection technique to the numerical simula- 
tion. In this figure, we have plotted the real and imaginary 
parts of the complex amplitude (|10l) for four g-modes of 
degrees £ — [1, 2] and radial orders n = [0, 1]. When a 
standing gravity wave occurs in a hydrodynamical simu- 
lation, its complex amplitude cg n behaves as 



ci n oc exp(-at) exp(iw^), 



(11) 



case of an excitation by an oscillating bubble dealt with in 
Paper I, the temporal Fourier transform of the mode am- 
plitude leads to a single peak centered around the eigen- 
frequency uJi n , with a width that is proportional to the 
viscosity (see Fig. 6 in Paper I). In other words, each g- 
mode obeys in this situation the same law as that of a 
linearly damped free harmonic oscillator. 

However, in the case of a random excitation by pen- 
etrating plumes shown in Fig. the real or imaginary 
parts now evolve either chaotically around zero when the 
mode is not excited, or in a periodic fashion as cos u>i n t 
or sinujint when the mode is excited. Indeed, some wave 
events are well visible, particularly for times t — [400—600] 
in the £ = 1, n = diagram, but the time evolution 
is mainly chaotic, suggesting that such wave events are 
difficult to extract. As a consequence, taking a temporal 
Fourier transform over the whole simulation of the mode 
amplitude C£ n (t) makes no sense as we will mix together 
wave and non-wave events. This is illustrated in Fig. ^5] 
where we computed the temporal Fourier transform of the 
mode amplitude in Fig.Elfor £ — 1 and n — 0. Comparing 
to Paper I, single peaks centered around the theoreti- 
cal eigenfrequencies toi n have been replaced by a forest 
of peaks roughly centered around wio, i.e. the mixing of 
wave events with non-wave events degrades the quality of 
the anelastic projection. Nevertheless, we remark that the 
spectrum in Fig. 1101 is better than the one obtained in the 
right-hand panel of Fig. [|J] with the classical method as 
peaks are now concentrated around the theoretical eigen- 
frequency. However, such spectra do not permit a detailed 
study of the amplitudes of g-modes propagating in the 
RZ, as other contributions (e.g., convective velocities in 
the overshoot region as well as penetrations of downward 
plumes) interfere with the time evolution of each mode 
amplitude. 

4.3. Time-frequency diagrams to extract wave events 

In order to accurately extract the hidden wave events, 
we use time-frequency diagrams of the mode amplitude 
cen(t), that is, the temporal Fourier transforms are com- 
puted by using a sliding window of fixed width (e.g., 
Flandrin & Stockier 1999). Assuming that this width is 
At (it is moreover beneficial to choose a multiple of the 
mode period) , we perform the following Fourier transform 
at time t 



c ln {t,L0) = I c ln (t') e iuJt ' dt', 



where uj£ n is the mode eigenfrequency and a a coefficient 
proportional to the diffusion process. For instance in the 



t+Ai/2 
t-At/2 

and thus iterate the process at the next time t + St, 5t 
being the timestep of the simulation. We then obtain a 2-D 
representation of the power spectrum |cf n (w)| 2 in a time- 
frequency plane (i, to) which highlights the time intervals 
during which the corresponding g-mode is really excited 
in the RZ. 

To illustrate the utility of this method, we focus on 
the g-mode at I = 1 and n — before applying it to other 
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Fig. 11. Time- frequency diagram of the amplitude of the 
g-mode at £ = 1 and n = using a temporal window of 
width At = 2Tio, with T w = 2tt/uj w ~ 22.6 the mode pe- 
riod. The horizontal dotted line corresponds to the mode 
frequency uiiq. 

modes; see Fig. 1111 Time intervals during which this mode 
is excited to significant amplitudes are well extracted, as 
many bumps appear along the line a->io, especially in the 
range t — [400 — 600] where three bumps are present. 
In order to isolate precisely and automatically the most 
powerful wave events, we apply the following procedure, 
illustrated in Fig. El still with our test mode at £ = 1 and 
n = 0: 

— we first compute a mean profile of the time-frequency 
diagram around the mode frequency ujg n (shown in the 
upper panel). 

— we then build what we call the "event function £" , that 
is, a function which is non-zero only when the previous 
mean profile is higher than its mean value: 



if /(*) < (/>* then 8 = 0, 
if /(*) > (f)t then 8 = 1, 



(12) 



and iterate the process four times by restarting from 
the new amplitude profile £ x /. We then obtain the 
final event function £ under a binary form (i.e. a suc- 
cession of and 1), shown in the middle panel in Fig. 
1121 Six events are clearly isolated, four of them being 
clustered in the range t — [250 — 620]. 
— finally, we apply the event function £ to the 
time-dependent amplitude C£ n (t). It thus emphasizes 
the time intervals during which cg n (t) behaves as 
exp(iu^ n i), that is, during which the corresponding g- 
mode is excited in the RZ; see the filtered real and 
imaginary parts of c\q in the lower panel of Fig. 1121 

We obtain the longest wave event at I = 1 and n = 
during times t — [344 — 435], which approximately corre- 
sponds to four mode periods, i.e. St = 91 ~ 4Tiq. This 
is confirmed in the snapshot of the velocity field at time 
t = 405 where a large-scale velocity field, signature of the 
propagation of the g-mode at £ = 1 and n — 0, is present in 
the bottom radiative zone (Fig. I13JI . Now, why do modes 
die out? The answer lies in Fig. El where we plot four 
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Fig. 12. Filtering of the mode £ = 1, n = 0: zoom in Fig. 
1111 around the frequency u>iq (upper panel); the resulting 
event function £ given by (|12|) (middle panel); the corre- 
sponding real (solid line) and imaginary (dotted line) parts 
of the mode amplitude after the filtering (lower panel). 



£=405 



o 

x 



Fig. 13. Snapshot of the velocity field superimposed on 
the internal energy fluctuations for time t = 405. Note 
the large-scale velocity field in the bottom radiative zone 
mainly due to the standing g-mode at £ = 1 and n = 0. 

snapshots of the velocity field for times t = [618 — 624] 
which correspond to the end of the events group clustered 
in the range t = [250 — 620]. The mode pattern is simply 
destroyed by the penetration of a plume which is born in 
the upper part of the CZ at time t ~ 615, then crosses 
the CZ and enters the RZ at time t = 620 where it kills 
the mode propagation. As a consequence, the large scale 
velocity field associated with the mode disappears and the 
event function £ becomes zero. 

Before generalizing this formalism to more modes, it 
is instructive to re-apply during this longest wave event 
the classical method discussed in H3.ll in order to com- 
pare both the spectrum and the simulated vertical mass 
flux with the theoretical ones. This is what we did in Fig. 
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penetrating plume entering in the RZ at time t ~ 620. 



1151 which is the same as Fig. at £ = 1 except that we 
focus now on times t = [300 — 450]. As expected, the 
bump around loio is more pronounced in the RZ (upper 
panel) and a comparison between the shark fin vertical 
profile deduced from a zoom around lu±o with the one 
computed from the theoretical anelastic mode shows an 
almost perfect agreement there (lower panel). It means 
that the dynamics of this region is well reproduced by our 
linear anelastic model. On the contrary, as g-modes cannot 
propagate in an unstable layer, it is normal to find a large 
discrepancy between these two profiles in the convection 
zone. 

We then apply this method to the g-modes of the first 
three degrees £ = [1, 2, 3] and radial orders n = [0, 1, 2], 
with resulting event functions £ given in Fig. 1161 That al- 
lows us first to show that the second bump located around 
u> ~ 0.2 in the upper panel in Fig. 1151 is due to the prop- 
agation of the mode £ = 1, n = 1, as the correspond- 
ing event function is not zero for times t ~ 400. Second, 
the assembly of the whole event function permits us to 
perform a statistical study of the mode lifetimes, that is, 
to compute the PDF of the duration of events (Fig. I17|) . 
This PDF is peaked around 2, meaning that the lifetime 
of a mode is approximately twice its period. Compared 
to the solar acoustic modes for which T osc ~ 5min and 
21if e ~ hour ~ 20T OSC , such a ratio is very small, i.e. 
g-modes patterns are rapidly destroyed by the following 
penetrating plumes and it may be a problem for their de- 
tection at the star surface. 

4.4. Kinetic energy due to g-modes 

Using the previous time-frequency diagrams for every g- 
mode allows us to find the time intervals during which 
IGWs propagate in the RZ. Now we want to quantify the 
efficiency of this stochastic excitation by using, say, some 
wave flux in the vertical direction. However, as we impose 



Fig. 15. Upper panel: same as in Fig. 0for £ = 1 except 
that the temporal Fourier transform has been computed 
for time t = [300 — 450] only. Lower panel: comparison 
between the corresponding vertical mass flux pw (solid 
line) and the theoretical one (dashed line) computed from 
the anelastic eigenvector at £ = 1 and n = 0. 

the reflective boundary condition w — for the vertical 
velocity both at the surface z = Z\ and the bottom z = 
z\ of our computational domain, we have standing waves 
rather than propagating waves and no flux is carried by 
waves in the vertical direction. As in Paper I, we thus chose 
to focus on the kinetic energy embedded in g-modes using 
the following Parseval type relation valid in the anelastic 
subspace 

/ pu 2 (x,z) dx dz = V |q„| 2 + "rest", (13) 

where the left-hand side corresponds to the total kinetic 
energy in the simulation at a given time t and the term 
"rest" contains all the contributions which are not due 
to IGWs (the demonstration is given in Paper I). This 
relation is very useful as it allows to quantify the kinetic 
energy embedded in every g-mode (£, n). Indeed, using 
the classical Parseval relation just allows one to quantify 
the amount of (say) kinetic energy contained in a mode 
of given degree £ and contributions coming from g-modes 
as well as p-modes or any turbulent motion are mixed 
together. The advantage of our anelastic subspace relation 
(|13|l is that it lifts this degeneracy in the radial order n by 
isolating the g-mode contributions. Of course, this relation 
should be applied only during times when IGWs effectively 
propagate, i.e. when a n cx exp(iw£ n t) or £ = 1, such that 
the \ci n \ 2 contributions make sense. 

The temporal evolution of the total kinetic energy 
-Etot embedded in the simulation (i.e. the left-hand side in 
Eq. IT3)l is illustrated as the solid line in Fig. IT%k . while the 
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Fig. 16. The event functions £ for g-modes of degrees £ — 
[1, 2, 3] and radial orders n = [0, 1, 2], with their associated 
periods Tt n . 
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Fig. 18. a): time evolution of the total kinetic energy 
(solid line) embedded in the simulation, with its compo- 
nent which is only due to g-modes (dot-dashed line), b): 
ratio between the two. 

t ~ [400 — 600] when the 1 = 1, n = mode is strongly 
excited. This large ratio is interesting, since it means that 
internal gravity waves may contain a non-negligible part 
of the total kinetic energy of the coupled system formed by 
the two neighboring convective and radiative zones. This 
is a direct demonstration that the excitation of IGWs by 
penetration plumes can be quite efficient, at least in 2-D. 



0.8 




lifetime of wave events 

Fig. 17. PDF of the lifetime of wave events, in units of 
the mode periods. 

dot-dashed line in the same panel corresponds to the con- 
tribution £"igw coming from g-modes only (the right-hand 
side sum in Eq. ^3 where the I = [1 — 3] and n = [0 — 2] 
modes have been considered). The interesting quantity is 
of course the ratio between the two, that is Ejqw / E tot , 
which emphasizes the efficiency of the IGW excitation 
by the downward plumes (Fig. 118b). It emerges that g- 
modes contribute up to about 40% of the total kinetic 
energy when they are excited, for example in the range 



5. Conclusions 

We have investigated numerically the excitation of inter- 
nal gravity waves by the penetration of convective plumes 
into an adjacent stably stratified zone. This problem is 
intimately related to the transport processes of chemicals 
and/or angular momentum in radiative zones of cool stars, 
such as the Sun. The knowledge of both spectrum and am- 
plitudes of the internal waves field is crucial. 

After recalling the main features of 2-D hydrodynam- 
ical simulations of penetrative convection, we focused on 
the problem of the mode identification by first using a 
classical method based on successive Fourier transforms 
of the vertical mass flux over the whole simulation. We 
thus showed that this tool is not adapted to detect the 
g-modes propagating in these simulations as the resulting 
spectra are very noisy, preventing a quantitative study of 
the phenomenon. 

We then introduced our method for detecting accu- 
rately the internal waves in the radiative zone. It is based 
on the projection of the simulated velocity field onto the 
anelastic eigenmodes that are solutions of the associated 
linear eigenvalue problem for the perturbations. Indeed, 
when a standing wave is present near the bottom radia- 
tive zone, its spatial structure is that of an eigenmode and 
the coefficient of the projection onto this eigenmode gives 
the mode amplitude. This amplitude is of course time- 
dependent as the internal wave is only generated after the 



Boris Dintrans et al.: IGW excitation by penetrative convection 



11 



penetration of plumes, and is then dissipated under the ac- 
tion of both viscous and radiative diffusions. This leads us 
to use time-frequency diagrams to isolate the most pow- 
erful wave events and to construct what we called "the 
event function £", that is, a binary function (0/1) which 
is set to 1 only when the mode amplitude is higher than a 
given threshold. As an example, we focused on the g-mode 
at I = 1 and n — 0, whose period is Tio = 22.6. We ex- 
tracted six wave events in our particular 2-D simulation, 
the longest one corresponding to four mode periods. We 
then showed the intricate link existing between the mode 
and the downward plumes as they can either excite it or 
destroy it! 

The extension of this study to the g-modes at t = [1—3] 
and n = [0 — 2] allowed us to compute the PDF of the 
mode lifetimes. We found that the mean lifetime is only 
around twice the period of the mode. The shortness of this 
lifetime may be a problem from an observational point 
of view where one needs lifetimes as big as possible (the 
large-scale solar acoustic modes have lifetimes of about 
the day, i.e. several hundreds of times their mean period). 

Finally, we looked at the kinetic energy content of the 
excited g-modes and showed that up to 40% of the total 
kinetic energy at times may reside in g-modes. This level 
is only reached during a fraction of the time, and the mode 
kinetic energy varies considerable with time. Nevertheless, 
when the modes are excited, the corresponding velocity 
field in the radiative zone has an amplitude that may lead 
to an efficient wave transport there (through the advective 
term M wavo -V). 

It is clear that our detection method allows a quan- 
titative analysis of the problem of g-mode excitation by 
penetrative convection. Following this work, we have been 
doing a parametric study of the influence of the convective 
flux on the mode amplitudes, by trying to predict these 
amplitudes from mixing-length arguments. Some recent 2- 
D simulations by Kiraga et al. (2003) indeed suggest that 
such mixing length models systematically underestimate 
the strength of the internal wave field. 

In the same way, the depth of the penetration in the 
stably stratified zone is without doubt a key parameter 
in the excitation mechanism by penetrative plumes. This 
penetration strongly depends on the local value of the 
Peclet number Pe: (i) large values of Pe mean that the 
advection is greater than the radiative diffusion such that 
the plume keeps its identity and is stopped very rapidly 
by the buoyancy braking, leading to a tiny penetration; 
(ii) small values of Pe mean that the plume thermalizes 
very rapidly with its surrounding and the buoyancy brak- 
ing disappears, leading to a large penetration. It would 
then be interesting to further study the influence of the 
Peclet number on the mode amplitudes, by computing for 
example a grid of 2-D polytropic models with different in- 
dexes 7713, in order to play with the Peclet number jump at 
the base of the convection zone (Eq. EJ. Likewise, the in- 
fluence of the dimensionality; i.e. the differences between 
2-D and 3-D also need be investigated. 
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Appendix A: The initial stratification in the three 
layers 

The initial vertical stratification is computed using polytropes 
of various indexes for which 

Pl + l/m mm 
oc p ' or p rxT , 

where m denotes the polytropic index. Such polytropic solu- 
tions are useful in numerical simulations of convection as they 
allow to easily specify the stability (or not) of a layer to con- 
vection. Indeed, the well known Schwarzschild criteria for con- 
vection are (e.g., Hansen & Kawaler 1994) 

f V < V ad => STABLE, 

{ (A.l) 
I V > V ad => UNSTABLE, 

where V = d\nT/d\nP and V ad = 1 — I/7 is its value in 
the case of an adiabatic stratification. The using of polytropic 
solutions leads to the following simple form for the "del" 



meaning that a polytropic layer of index m will be stable or 
unstable to convection following 

f V < V ad => m > m ad : STABLE, 

{ (A.2) 
I V > V ad => m < m ad : UNSTABLE. 

Since 7 = 5/3, we have V ad = 2/5 and m ad = 3/2 in this case 
and a polytropic layer will be convectively stable if m > 3/2 
and unstable if m < 3/2. 

Once the three polytropic indexes of the superposed layers 
have been fixed we first compute the corresponding radiative 
conductivity profile by assuming that all of the energy flux 
.Fbot that we impose at the bottom is initially transported by 
radiation, that is, 

£i = ^%-l)(l + 0> (A-3) 
9 

in the layer i (i = [1, 2, 3]). In this formalism, the radiative con- 
ductivity is thus a constant in each layer and its three different 
values are joined using polynomials of the third-order. 

In Fig. lA.lfa we show an example of such a profile for a 
polytropic layered system with indexes mi = —0.9, m% = —0.5 
and ms = 3. As expected, the radiative conductivity in the CZ 
is below the critical value /C ad given by 

£ad = ^(7-l)(l+mad), 

that is, the value deduced from Eq. 1A.3I using m = m ad . More 
surprising could be the value /Ci < /C ad in the surface layer as 
one rather expects a stably stratified layer here, as observed in 
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Fig. A.l. Initial vertical profiles for a layered system with 
m = [—0.9,-0.5,3] of the radiative conductivity (a) and 
entropy (b). The vertical dotted lines delimit the convec- 
tion zone, whereas the horizontal dashed line in the (a)- 
plot corresponds to the critical value /C a d below which the 
layer is convectively unstable. The dot-dashed line in the 
(b)-plot corresponds to the entropy profile in the relaxed 
state. 

solar-like stars where a very thin superficial stable layer exists. 
In fact, this layer is stable due to its efficient cooling [Eq. Q] 
such that the Schwarzschild criterion does not apply. We thus 
adopt a very weak value for /Ci, by taking mi = —0.9 in all 
of the simulations, and compute the initial hydrostatic strati- 
fication of this layer by assuming that it is already isothermal 
with de/dz = 0. 

Concerning the initial stratification of the CZ, things are 
more involved. Assuming an initial polytropic stratification 
with V2 = 1/(1 + m-i) is clearly not a good idea, as an effi- 
cient convection is always associated with an almost adiabatic 
stratification with V2 — V a( j. In this case, the relaxation of 
the CZ towards its adiabatic state would take a lot of numer- 
ical timesteps. One simple solution of this problem consists 
in starting from an adiabatic stratification in the CZ by im- 
posing V2 = V ac j. However, this solution does not take into 
account the entropy jump which appears at the top of the CZ 
(the difference between constant entropy in the CZ and the en- 
tropy minimum in the photosphere; see e.g., Abbett et al. 1997; 
Ludwig, Freytag & Steffen 1999) such that the relaxation time 
would stay important. 

One solution of this relaxation time problem consists in 
starting from a mixing length stratification where the (local) 
superadiabatic gradient in the CZ is modeled using the follow- 
ing mixing-length argument 



v 2 = v ad + a 



pel 



2/3 



where c s — \J 7(7 — l)e denotes the local sound speed, a ~ 1.5 
is a free mixing-length parameter and -F conv is the convective 
flux given by 



Fccmv — Fi 



bot 



^rad — ^1 



bot 



l--( 7 -l)(l + m 2 ; 
9 



Finally, the initial stratification of the bottom RZ is sim- 
ply computed by assuming that all of the bottom flux is 
transported in this layer by radiation and one sets V!{ ad = 
1/(1 +ma). 

To summarize, the initial stratification of our three-layer 
model is computed from 



dlnp 

dz 



7- 



1 e 1 



where e obeys to the following set of equations in the three 
layers 







de 
dz 

de _ g 

dz 7 — 1 



if Z\ < z < Z-2, 

Vmlt -r „ 
2 it Z2 < Z < 23, 



de g ra d . e , , 

— = -V 3 if Z 3 < Z < 24. 

K dz 7 — 1 



and these differential equations are iterated until the density 
at the BCZ is equal to 1; i.e. we impose p(z — Z3) = 1. 

In Fig. O ) we show the resulting vertical profile of 
the initial entropy s/c p = ln(p/p 7 )/7 for the same indexes 
m = [—0.9, —0.5, 3.]. The solid line denotes the initial entropy 
whereas the dot-dashed line corresponds to its profile in the 
relaxed state. First, as the entropy gradient in an isothermal 
layer is a constant, one verifies that s tx z for z = \z\,zi\- 
Second, the comparison between the initial and relaxed profiles 
shows that the mixing-length stratification well reproduces the 
entropy jump at the top of the CZ: the strong mixing taking 
place in the deep layers of the CZ leads to an almost flat en- 
tropy profile, which disappears at the base of the photosphere. 
As a consequence, the computing time needed to relax towards 
this solution is considerably reduced by using mixing-length 
solutions. 
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